% Script Eg11_1
% Examines the accuracy of the simple daylight-predictor model.
% Assumes that InsightData is a subdirectory of the working directory.

% Data is available for these cities..
Cities =  {'Anaheim',  'Anchorage', 'Arlington', 'Athens', 'Atlanta',...
    'Baltimore','Bangkok', 'Beijing', 'Berlin', 'Bogata', ...
    'Boston', 'BuenosAires', 'Cairo', 'Chicago', 'Cincinnati', ...
    'Cleveland', 'Denver', 'Detroit',   'Honolulu',  'Houston', ...
    'Ithaca', 'Johannesburg', 'KansasCity', 'Lagos', 'London', ...
    'LosAngeles', 'MexicoCity', 'Miami', 'Milwaukee', 'Minneapolis',...
    'Moscow',     'NewDelhi',     'NewYork', 'Oakland', 'Paris',...
    'Philadelphia', 'Phoenix', 'Pittsburgh', 'RiodeJaneiro', 'Rome',...
    'SanFrancisco', 'Seattle', 'Seoul', 'Sydney', 'Tampa',...
    'Teheran',  'Tokyo',  'Toronto', 'Washington', 'Wellington'};

% Choose a subset of the city indices 1 through 50 to investigate...
CityChoices = [4 21 25 33 50];
for k = CityChoices
    % The sunrise/sunset data for city k...
    fn = [ pwd '\InsightData\11_1\' Cities{k}  '.dat'];
    % Get the actual and predicted daylight values and display...
    [Name,Lat,Long,Hours] = GetRiseSetData(fn);
    HoursEst = Daylight(Lat);
    err = ceil(max(Hours-HoursEst)*60);
    Days = (1:365)';
    %    Display...
    plot(Days,Hours,'k',Days,HoursEst,'r')
    legend('Actual','Model')
    title([Name sprintf('  Latitude = %4.1f   Max Error = %1d minutes',...
                        Lat,err)])
    ylabel('Hours of Sunlight')
    axis([0 365 4 20 ]);
    set(gca,'Xtick',[31 59 90 120 151 181 212 243 273 304 334],...
            'XTickLabel',[]);
    grid on
    z = 3.5;
    text(7,z,'Jan');   text(37,z,'Feb');  text(67,z,'Mar')
    text(97,z,'Apr');  text(127,z,'May'); text(157,z,'Jun')
    text(187,z,'Jul'); text(217,z,'Aug'); text(253,z,'Sep')
    text(283,z,'Oct'); text(313,z,'Nov'); text(343,z,'Dec')
    shg
    pause(1)
end